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Abstract 

Sampling  techniques  and  exact  solutions  of  Riemann 
Problems  are  us^  in  a random  choice  method.  This 
procedure  is  used  to  obtain  the  numerical  solutions  of  a 
q.stem  of  conservation  laws  which  describes  the  dynamics 
of  flow  for  small  amplitude  twx>-dimcnsional  shoclcHiaves. 
An  intrinsic  coordinate  system  is  used  to  formulate  the 
model 

1 Introduction 

Aamracy  of  numerical  solutions  and  effidcncy  of 
numerical  sdiemcs  arc  major  concerns  in  obtaining 
numerical  solutions.  Moreover,  the  numerical  solution 
a:  the  jump  discontinuities  called  shocks  should  remain 
sharp,  stable  and  traitsports  the  discontinuities  at  the 
correct  physical  speed.  Random  variables  have  been  used 
to  control  numerical  dissipation  or  to  control  numerical 
viscosity.  Basicalfy,  random  variables  appear  either  as  a 
component  added  to  the  deterministic  equation  to  study 
the  effect  of  numerical  viscosity  or  they  arc  used  to 
sample  the  solution  at  a randomly  chosen  point  to  obtain 
a numerical  solution  which  preserves  some  mathematical 
propenies  of  the  solution  function.  The  purpose  of  this 
paper  is  to  present  a random  choice  method  for 
computing  the  numerical  solution  of  twxr-dimcnsional 
small  amplitude  sbockwav^cs.  The  numerical  random 
sampling  procedure  is  a shock  capturing  and  a marching 
time  method  for  solving  ^stem  of  conservation  laws. 
The  random  sampling  procedure  consists  of 
approritnating  the  numerical  solution  by  a piecewise 
constant  state  at  each  time  step  and  proceeding  to  the 
next  time  step  by  sohing  the  corresponding  problems 
formed  by  the  constant  on  the  neighboring  spatial 
intervals.  It  is  well-known  that  the  exaa  solution  of  non- 
linear system  of  partial  differential  equations  arising  in 
fluid  flow  problems  even  with  smooth  initial  data 
develops  shocks  (jump  discontinuities;  in  a finite  time 
intetvtd.  Thus  it  is  not  unnatural  to  approrimate  their 
initial  data  with  consuni  states. 

The  sampling  procedure  is  based  on  approximating  the 
numerical  solution  of  the  gnen  problem  with  a sequence 
of  elemcntaty  problems,  known  as  the  Riemann 
problems.  These  Riemann  problems  can  be  thought  of 
as  information  source  about  the  solution  within  each 


spatial  mesh  interval  More  importantly  thigr  provide 
valuable  information  on  wave  interaaion. 

Godunov  jl]  initiated  utilizing  the  solutions  of  the 
Riemann  problems  as  building  blocks  for  the 
construction  of  numerical  solution  of  the  nonlinear 
hyperboUc  partial  differential  equations.  Godunov 
replaced  the  initial  data  by  a piecewise  constant  states 
with  jump  discontinuities  at  the  middle  of  spatial  mesh 
interval.  Then  the  exact  solution  of  this  Riemann 
problem  at  the  first  time  step  is  calculated.  To  proceed 
to  the  next  time  step  replace  this  exaa  solution  a new 
piecewise  constant  state  approximation  and  soho  the 
corresponding  Riemann  problem  and  maintain  integral 
properties  of  the  conserve  variable. 

Another  utilization  of  Riemann  problems  in  obtaining 
the  solution  of  conservation  laws  was  initiated  1^’  Glimm 
(2]  who  followed  Godunov  as  far  as  obtaining  the  exaa 
solution  of  Riemann  problem  and  then  the  value  of  the 
new  approximated  solution  at  the  new  time  step  is  taken 
to  be  the  exaa  solution  evaluated  at  a random  point  on 
that  mesh  interval  Ibis  solution  is  conservative  on  the 
average,  however,  has  the  advantage  that  near  jump  the 
solution  is  inacmcnied  dther  by  the  amount  of  jump  or 
not  at  ail  This  forces  that  an  initially  sharp 
discontinuities  remains  sharp.  Chorin  (3]  developed 
Glimm’s  random  choice  method  into  a numerical 
technique  The  random  choice  method  by  its  way  of 
construaion  propagates  shocks  without  introducing  any 
dissipation  and  the  method  is  unconditionally  stable 
However,  because  of  approrimating  solution  at  a 
randomly  chosen  point  a small  amount  of  statistical  noise 
enters  into  the  solution  which  is  acceptable  within  the 
accuracy  imposed  by  discretization  of  model  problem. 

2 Two-Dimensional  Flow  Problem 
The  equations  describing  the  two-dimensional  flow  of 
shockwaves  with  a source  term  in  fluid  dynamics  for 
compressible  fluid  may  be  written  in  the  form 

(1)  Uj.  + flu)^+g(u)y  = b{u,Xsy,t) 

where  f and  g arc  physical  fluxes,  h is  a source  term  and 
the  unknown  quantity,  u is  a funaion  of  x,  y,  t.  Denoting 
the  front  coordinate  by  a and  letting  the  coordinate  B be 
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the  arc  length  measured  from  a reference  point  along  the 
front,  then  the  successive  front  positions  are  given  by  the 
family  of  curves,  a = constant  and  the  ray  positions  by 
the  family  of  curves,  B = constanL  By  using  this  intrinsic 
coordinate  system  a fi  (see  Whitham  {8])  where  a and  6 
are  functions  of  x,  y,  t,  equation  (1)  can  be  written  as 

(2)  w,  + F{w)p  = G(w,  a,  P) 
subject  to  the  initial  condition  ^en  by 

(3)  w-(0,  p)  =w;j(P)  . 

Equations  relating  x and  y to  a and  B are  given  by 

= (1  +Viin<J))  cos  (0)  = -Asin(0) 


y,  = (1 +1/4m4>)  sin(0)  yp=Acos(0) 

Here  8 is  the  angle  that  each  front  makes  with  the 
positive  x-axis,  A is  the  cross-sectional  ray-tube  area,  and 
m is  the  acoustic  Mach  number.  For  small  amplitude 
twx>-dimcnsional  shockwaves  we  have 


i 

f A ^ 

f ] 

AQ 

F(w)  = 

(m(J)-02)/2 

. 0 , 

(->) 


G(w) 


^ 0 ^ 

0 

r 4CZ, 


where  C is  the  local  sound  speed,  ^ is  the  nonlinearity 
constant  which  depends  on  the  media  and  Z is  the  area 
under  the  initial  pulse.  For  a detail  discussion  of  these 
equations  see  Zakeri  [4]-[5).  To  solve  (2)-(3)  we  use 
operator  splitting  method  to  remove  the  inhomogeneous 
term  G(w.  a,  B).  That  is,  first  we  sohe  the 
corresponding  one-dimensional  homogenous  problem. 


(5)  w,  + F(w)p  = 0 

by  sampling  procedure  and  then  we  use  its  solution  to 
determine  the  value  of  the  inhomogeneous  terra,  G(w,  o, 
B).  Finally,  we  solve  the  corresponding  ordinaiy 
differential  equations  (ODEs)  given  by 

(6)  = G{w,  a.  P)  . 


To  sohe  (6)  we  use  a common  ODE  sober  such  as 
Rungc-Kutta  or  a multi-level  method. 


g^n  by  (2)-(3).  One  of  the  advantages  of  formulation 
of  a model  problem  using  geometric  shock  dynamics  is  its 
simplicity.  To  develop  a random  choice  method  first  we 
must  define  a random  variable  defined  over  closed 
interval  V^J.  It  is  absolutely  necessary  that  the 
successive  values  of  the  random  variable  tend  to 
approximate  e«^ui-partitioning  the  closed  interval  [ 1*^) 

(see  Glimm  [?]).  To  generate  such  random  variable  let 
us  consider  a sequence  of  pseudorandom  integers 
generated 


0)  (iTDdic) 

where  k is  an  odd  posith’e  integer  and  Nq  is  an  arbitrary 
integer  less  than  k.  Let  us  define  an  equidistribuied 
sequence  random  variables,  o„  on  the  interval  [-V6  , 16] 
^ en  by 

1 

(8)  o„  = -2  - — . 

We  introduce  front-ray  grid  defined  by  mesh 

lengths  aa  and  aP  . The  solution  of  (2)-(3)  is  to  be 
calculated  both  at  grid  points,  Le.,  at 

P(n.  j)  = (naa  , jaP) 
and  at  the  center  of  rectangle  grid  point,  Le.,  at 

P{n+'A,  j+'A)  = ((j3+%)aa,  (j-»-%)aP) 

where  n and  j are  integers.  We  denote  the  approxirnate 


value  ofw  at  the  grid  point  by  w?  = (riaa  , jaP)  . 

Following  the  outline  given  above,  let  us  consider  the 
corresponding  local  Rictnann  problem  to  (2)  when  from 

is  at  riaa  along  with  the  piecewise  constant  initial  data 
given  by 

(9)  F,  + F(R)  p = 0 


F(riaa, P) 


; P < Jap 
; p i JaP 


where 


J-  = j-%(-D* 


c = + sort  (o^^i)  -1) 

Le.  e = 0 or  1 whenever  is  negathe  or  non- 
negative  respectively.  The  Rjeraann  problem  here  is 

sampled  at  (j+%)aP  and  at  (j-\^)aP. 


3 Numerical  Scheme 

We  develop  a numcricai  scheme  to  compute  the 
successive  shock  fronts  using  geometrical  shock  dynamics 
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When  is  non-negath-e  the  initial  data  for  Riemann 
problem  formed  by  using  information  at  grid  points 
P(n,j)  and  P(n,j+1)  and  if  is  negative  then  the 
initial  data  is  constructed  by  using  information  at  grid 
points  P(n,j)  and  P(n,j-1).  At  point  P(n+V^  j+V5)  we 
define 

On  each  mesh  interval  we  get  a local  Riemann  problem. 
In  order  to  assure  that  the  waves  produced  by  this 
sequence  of  local  Riemann  problems  do  not  interact  we 
must  have 

(10)  -^C(l+1/S>in4>)  < 1. 

Aft 

This  important  requirement  is  known  as  Cbuiant- 
Friedrichs-Lewy  (C^)  condition.  If  inequality  (10) 
holds  then  we  can  combine  the  solutions  of  the  Riemann 
problems  (9)  into  a single  exact  solutiort 


4 SolutKHi  of  Riemann  Problem 
The  main  part  of  a random  choice  algorithm  is  obtaining 
the  solutions  of  a sequence  of  local  Riemann  problems 
cflldcntly.  The  solution  of  a Riemann  problem  consists 
of  three  elementary  waves,  a backward  shock  wave  or 
rarefaction  on  left,  a slip  line,  and  a forward  shock  or  a 
rarefaction  on  right  A slip  line  is  a discontinuous 
solution  separating  two  constant  states  sudt  that  the 
angle  of  flows  remain  the  same  on  beth  sides  of  the 
disconi>nuity  line  while  Mach  number  is  arbitrary.  Slip 
lines  are  one  family  solution  between  the  backward  and 
forward  waves.  Le.,  between  rarefactions  and  shocks.  To 
solve  the  Riematm  problem  (9)  we  follow  Lax  |6j.  Let  us 
corrida’  the  following  initial  data  for  system  of  cquauons 
in  (9) 


(H) 


Rin^a, P) 


p < JaP 
p ^JaP 


where  subscripts  1 and  2 refer  to  values  of  w just  behind 
of  and  just  ahead  of  the  discontinuity  rcspcahely.  If 
these  two  values  are  equal  then  the  solution  of  (9)  is  a 
constant  state  and  its  value  is  equal  to  the  value  of  initial 
data.  However,  if  these  two  values  arc  different  then  the 
initial  Jump  discontinuity  will  propagates  in  the  form  of 
a center  expansion  wave  and/or  a shock  (Le.,  jump 
discontinuity  satisfies  the  entropy  condition.)  or  a contaa 
discontinuity.  In  order  that  solution  converges  to  a 
unique  weak  solution  of  (9),  it  roust  satisfies  the 
Rankine-Hugoniot  jump  condition  and  the  Oleinik 
entropy  condition.  At  the  shock,  let  us  define  the  values 
of  R(a,B)  just  behind  of  and  just  ahead  of  the  shod:  by 


= lira  R(a  , P)  JZj  = lim  J?(a  , P) 

p-r 

The  jump  conditions  for  ^stem  of  equations  (4)  are 
given  by 

dp  62 
da 

^ ^ da  2(^282-^161) 

^2!!^ 

The  entropy  condition  is  given  by 

FiR^)  -F[R)  F{R^)  -F{Rj) 

R^-R  R^-R^ 

for  any  R between  R2  and  R,.  The  entropy  satisfaction 
is  a major  concern  for  numerical  approximation  of 
solutions  of  nonlinear  fluid  flow  problems.  This  simply’ 
means  that  the  computed  solution  comerges  toward  the 
correct  physical  solution  as  the  mesh  sizes  of  intervals 
along  a and  B approach  to  zeros.  The  above  inequality 
can  be  written  as 

£(i?)  =F(J22)  +-H  (R-K2) 

satisfying  the  following  inequality 

(F(£)  -F[R))  (i?i-J^)  iO 

where  E(R)  defines  the  chord  connecting  left  and  right 
limiting  points  across  the  shock.  The  entropy  related  to 
the  first  component  of  F is  given  fay 

A2-A 

for  any  A between  A,  and  A^,  and  6 between  0,  and  02- 
Similar  inequalities  hold  for  other  components  of  F. 

4.1  Rarefaction  Waves 

Rarefaction  waves  arc  two  families  of  solutions  curves, 
forward  and  backward  waves.  In  this  section  we  compute 
the  simple  rarefaction  waves  of  system  of  equations  (9) 
which  can  be  reformulated  in  the  form 

U^*H{V)  Up  = 0 

where  H(U)  = H(A,  0,  m)  is  3 by’  3 matrix  given  by 
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dioice  method  developed  her<  with  those  solution^ 
obtained  using  the  method  of  chaiactcristics.  Cbnsidef 
the  initial  condition 


1 

(O  -1  0 \ 

'A 

o 

o 

6 

H = 

2A 

0 0 

1 

1 2A 

The  eigenvalues,  p and  thier  corresponding  eigenvectors, 
e of  the  matrix  H are 


= 0 , ± Q - 

^ 2A 


1 


For  the  svstem  of  equations  in  (9)  the  simple  rarelaaion 
waves  arc  the  continuous  solutions  of  (9)  of  the  form 


(13) 


U(a.  3) 


if  |<ii(i2i) 
if  -&=n(v) 

C 


R,  ijf 


where  v is  an  integral  cutvc  of  the  vector  field  of  the 
corresponding  eigcnvtxtor  connecting  the  two  constant 
states  such  that  the  corresponding  eigenvalue,  p is 
increasing  between  this  two  constant  states  fiom  left  to 
right.  Since  the  matrix  H has  three  distinct  eigenvalues, 
there  are  three  possible  rarefaaion  waves  through  any 
ghen  state.  These  rarefaaion  waves  are  the  integral 
curves  of  the  veaor  field  defined  by  each  cigcnveaor  of 
matrix  H,  Lc.,  each  cigcnveaor  is  tangent  at  each  point 
of  integral  curve.  Thus  for  the  eigenvcaor  e 
corresponding  to  the  eigenvalue  p of  matrix  H,  the 
integral  curves  arc  solutions  the  following  system  of 
equations 


dA  _ dd  _ dm 
1 -2\i^A/^ 


If  p = 0 then  the  integral  curves  of  its  corresponding 
cigcnveaor,  c = (1, 0, 0)  arc  curves  where  0 and  m arc 
both  constants.  Hence  a simple  rarefaaion  wave  of  the 
form  of  (13)  exists  if  the  left  and  right  values  of  0 and  m 
across  the  shodc  arc  equal,  in  addition,  p must  be  an 
inacasing  funaion  of  0 and  m from  left  to  right  across 
the  shock.  Therefore  there  is  no  A-rarefaaion  wave. 
0-raic&ak»  waves.  If  p is  not  zero  then  the  integral 
curves  of  its  corresponding  cigcnveaor  arc  curves  where 
m^A  is  constanL  There  arc  two  families  of  curves  where 
0 is  cither  positive  or  negative  along  each  integral  curve. 


5 Numerkal  Experiments 

We  compared  the  numerical  solutions  using  the  random 


A(0,  3)  =1 

tn(0,  3)  =0.01 


6(0,  3)  = 


^3  -!<3<1 

0 otherwise 


toother  with  equations  (2)  and  (4).  The  results  obtained 
with  both  methods  were  very  doK  to  each  other.  Th<7 
nuniericai  calculations  show  that  the  convergence  of 
solution  toward  the  exaa  solution  is  independent 
choice  of  the  odd  number  k in  (7)  as  long  as  k 
boutrded.  In  addition,  the  method  is  numerically  stab]^ 
All  the  striking  physical  features  of  ^stem  of  equation^ 
(2)  and  (4)  is  observed,  Lc.,  as  the  initial  frotri 
propagates  into  the  rest  state  the  center  of  hump 
becomes  flat  and  this  flat  re^on  propagates  on  both 
dircaions  untill  the  front  becomes  a flat  surface. 


6 Conclusion 

The  numerical  solutions  show  that  the  method  is  stable 
and  correaly  desenbes  the  important  physical  feature  of 
the  solution  of  the  model  problem.  The  various  dioic^s 
cf  random  number  generators  do  not  have  any  effea  oh 
the  accuracy  of  computed  solution  as  long  as  the  raodojh 
variable  tent  toward  the  equipartitionieg  of  the 
intervaL 
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ABSTRACT.  T infar  ccmpaitmentai  systems  have 
conqmtments  with  flows  to  and  from  the  compartments. 
Of  interest,  is  the  estimation  of  the  constants,  0, 
goverring  the  flows.  In  the  particular  sy^em  considered, 
only  one  c^n^nrlment,  out  of  several,  is  obsoved  for  n 
cases  over  k time  points.  A stochastic  model  is  used  with 
a maxicmm  likelihood  ^>proach  taken  to  the  estimation  of 
9.  The  algorithm  involves  iteratively  using  an  estimate  of 
9 to  solve  differential  equations  which  describe  the 
system,  and  inqirove  on  the  estimate  of  6 by  adding  a 
constant  multiple,  or.  of  an  increment  bi,.  .'  lien  s results 
rre  incorporated  to  obtain  required  deri\'atives.  Due  to 
non-zero  correlations,  a modification  to  Jetuurich  and 
Moore's  results  is  made,  involving  using  both  the 
observations  and  thei:  cross-products,  to  obtain  o6.  a is 
determined  with  Fletcher's  method.  A program  in  Turbo 
Pascal  inqilenrents  the  algorithm. 

INTRODUCTION.  Compartmenlal  systems  have  long 
been  a useful  tool  in  pharmacokinetics  O^^gner  (1971)). 
The  body  is  thought  of  as  a series  of  conqjartrocuts,  with 
a drug  moving  between  any  of  the  conqiartments.  For 
exanqile,  Gladtke  et  al.  (1979)  p.  36.  suppose  that  the 
body  may  be  n^resented  as  three  compartments  - plasma, 
muscle  and  extravascular.  .An  initial  muscle  injection  is 
given  and  the  levels  of  the  drug  in.the  plasma  are 
monitored  from  time  to  time.  The  drugs  will  flow  from 
nmscle  to  plasma.  .Additionally,  flow-  will  be  between  the 
extravascular  ^-stem  and  plasma,  and  from  plasma  to  the 
outside,  as  dqiicted  in  d.iagram  1 . 

Diagram  1: 
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traditional  deterministic  model,  these  are  the  rate 
parameters). 


Diagram  2 : 
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It  is  assumed  that  a bolus  injection  is  given  in  compartment 
1 and  only  compartment  2 is  observ-ed.  at  times  t,,t2,...,ti. 
Let  be  the  initial  concentration  in  compartment  I and 
C(t)  be  the  concentration  at  time  t.Then 
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STOCILASTIC  MODEL  . For  the  ’particle  model’,  as 
discussed  by  Purdue  (1974a),  it  is  assunred  that  there  are 
N particles  in  the  system  acting  indqrendently.  Tranritiems 
between  compartments  follow  a Markov  process,  with  the 
transition  probability  being  constant.  The  resulting  system 
of  equation  is  as  follows: 

dP(t)/dt  = P(t)A^  (I) 

where  P(l)  = (p^t))  is  a nonsingular  3x3  matrix  with  P5(l), 
ij=  1,2,3,  the  probability  of  a particle  transferring  from 
compartment  ’i’  to  compartment  ’j’  in  time  t and 


Outside! 


For  simplicity,  this  paper  will  be  restricted  to  a three 
COmpartmental  system  of  this  type,  with  conqiartments 
1,2,3  and  the  outside,  as  conqrartment  0;  flows  have  a 
parameter  attached  to  them  as  in  diagram  2,  (In  the 
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Under  these  assumptions  it  can  be  sho^^iii  that  the  exact 
distribution  of  C(t),  given  C(0)  is  a multinomial  (since  only 
C,(0)  is  nonzero),  but  the  distribution  of  C(t,),  given  Cft^, 
is  a convolution  of  multinomials.  Thus,  to  write  down  the 
exact  distribution  of  the  concentrations  in  compartment  2 for 
tj,t^...,t^  when  k is  large  is  impossible  since  it  w-culd 
consist  of  many  sums  whose  limits  are  complex.  Using  a 
diffusion  approximation,  it  has  been  demonstrated  that  the 
concentrations  over  time  have  a multivariate  normal 
distribution,  Lefaoczky  and  Gz^er  (1977).  Hence,  the 
distribution  of  the  concentrations  in  compartment  2 have  a 
(marginal)  multivariate  distribution,  with  mean  and 
variance-covariance  matrix  the  same  ^ that  given  by  the 
particle  model . Sicqrson  (1988)  has  snown  that  an  e^mator 
derived  from  the  maximum  likelihood  equations,  using  this 
^rpmximate  normal  -iisttibubon.  Is  still  a consistent 
asynqrtoticaliy  nomsl  ’>r,  if  the  particle  model  is  the 

correct  one.  Thus,  an  algcril.in  was  written  to  estimate  8 
using  the  approximsn;  aormai  disiribution. 


Suppose  that  is  the  concentration  in  compartment  2 for 
person  i,  i=l,2,...n  at  time  t,.  j = 1.2,...k.  3^  = 
(Xj„Xs,....Xa)  are  the  k observations  of  ith  case.  X„X-...., 
Xj  are  indqiendently  and  identically  distributed  as  a normal 
variable  with  mean  and  variance  which  are  functions  of  the 
probability  matrix  P.  By  considering  a vector  comprised  of 
the  sufficient  statistics  of  the  covariance  of  X,  it  era  be  seen 
that  a normally  distributed  random  variable  is  ft  licesr 
exponential  and  so  the  algorithm  of  Jennrich  and  Moore 
(1975)  ,hence  referred  to  as  JMA,  can  be  used  to  Ssd  the 
approximate  maximum  likelihood  estimator  . We  Kke  the 
k,xl  vector  Y,  where  k,=  k(k-r3)/2  and 
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THE  ALGORITHM.  Suppose  that  the  k,x5  matrix  of 
partial  derivatives  with  respect  to  i=0,1,..,4  is  denoted 
by  641/68.  According  to  the  JMA,  for  a given  0,  replace  $ 
by  •+a55, where 

This  requires  that 

(1)  P(.)  and  its  derivatives  with  re^ject  to  0 be  obtained 

(2)  06  be  calculated. 

(3)  An  appropriate  or  be  found. 


m P(-)  and  its  derivatives: 

For  this  conqiartmental  ^stem.  P(.)  and  its  denvatives  c^ 
be  calculated  analyJcally.  How-ev'er,  to  maintain  generality 
for  the  program,  it  was  decided  to  obtain  I^.)  and  its 
derivatives  by  applving  the  method  developed  fay  Allen 
(1987)  to  columns  of  P(.). 

The  steps,  for  each  column  of  r,  are  as  follows: 

(i)  Fiiid  R so  that 

- R^AR  = T , a triangular  mafrix. 

- R^  « I,  the  identity  cztrix 


(ii)  Using  R,  obtain  a triangular  system  of  equations; 

so,  a backward  living  technique  can  be  used,  with 
the  initial  condition  90)  - I .-satisfied. 


3r  any  a = 5 ,2.5.4,  liie  equation  (1)  becomes 

( t)  ^TV\  ( C)  *R  R t)  , 

at  * 00, 
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and  Allen's  method  can  be  tised  again  to  solve  for  V.(t), 
with  V,(0)  = 0,  as  the  initial  condition 


E^ric-a.E^ric-iA«* 


Using  properties  of  normality,Y  has  a mean,  41(8),  and 
variance-covariance  matrix,  E(fl),  which  is  a function  of  the 
mean  and  variance-covariance  matrix  of  X,  and  therefore  of 
P,  also. 


(2)  Calculation  of  bOz 

Using  Wilkinson’s  algorith.ms,  Wilkinson  (1971),  a 
Qioleskv*  decomposition  is  used  to  Bnd  where 
E^(2^)^=E.  Using  this,  the  following  equation  is  solved  to 
obtain  o9: 

l»i  = 


dd 
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(31  An  or. 

Fletcher's  descent  method.  Fletcher  (p.26,  1980)  is  used  to 
find  a.  It  uses  the  first  derix-atives  only.  Assuming  uniform 
continuity  conditions,  it  will  achie\'e,  at  least,  a local 
optimum  if  an  optimum  exists  and  if  the  starting  value  is 
close  enou^. 

SUMMARY.  A stochastic  approach  rather  than  the 
traditional  detenniiustic  model  approach  is  taken  to  a 
particular  conqiartmenial  noodel,  where  only  one 
con^sartment  is  observed.  An  algorithm  is  developed  which 
uses  the  JMA  to  obtain  maximum  likelihood  estimates  from 
a diffusion  ^proximation.  The  program  is  written  in  Turbo 
Pascal  and  can  be  generalised: 

(i)  to  other  linear  compartmental  s)'stems, 

(ii)  for  u to  be  fimctions  of  time. 

Work  is  being  done: 

(iii)  to  incorporate  measurement  error  in  the  nKxiel 

(iv)  to  include  people  variation 

An  important  stqr  for  its  general  use  wxruld  be  to 
incorporate  this  program  in  a general  pharmacokinetic 
program  so  that,  in  a user  friendly  enviromnent,  its 
estimates  could  be  easily  compared  to  those  obtained  from 
other  nrethods. 
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